Usenet 1993 July/InfoMagic USENET CD-ROM July 1993.ISO/sources/unix/volume10/complex-lib

Text File | 1987-07-14 | 24.2 KB | 1,106 lines

Path: uunet!rs
From: [e-mail] (Rich Salz)
Newsgroups: comp.sources.unix
Subject: v10i056:  Complex arithmetic library
Message-ID: <[e-mail]>
Date: 16 Jul 87 00:12:15 GMT
Organization: UUNET Communications Services, Arlington, VA
Lines: 1095
Approved: [e-mail]

Submitted-by: [e-mail] (Doug Gwyn)
Posting-Number: Volume 10, Issue 56
Archive-name: complex-lib

[  Here is a package to do complex arithmetic; blame me for the Makefile.
    --r$  ]

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of shell archive."
# Contents:  Makefile complex.3 complex.h cx_test.c cxadd.c cxampl.c
#   cxconj.c cxcons.c cxcopy.c cxdiv.c cxmul.c cxphas.c cxphsr.c
#   cxscal.c cxsqrt.c cxsub.c
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f Makefile -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"Makefile\"
else
echo shar: Extracting \"Makefile\" \(395 characters\)
sed "s/^X//" >Makefile <<'END_OF_Makefile'
XALL=complex.3 complex.h libcomplex.a
XOBJS=\
X    cxadd.o cxampl.o cxconj.o cxcons.o cxcopy.o cxdiv.o cxmul.o \
X    cxphas.o cxphsr.o cxscal.o cxsqrt.o cxsub.o
X
Xall:    $(ALL)
X
Xinstall:    $(ALL)
X    @echo install $(ALL) according to local convention.
X
Xcx_test:    cx_test.c libcomplex.a
X    $(CC) $(CFLAGS) -o cx_test cx_test.c libcomplex.a
X
Xlibcomplex.a:    $(OBJS)
X    ar r libcomplex.a $(OBJS)
X
X$(OBJS):    complex.h
END_OF_Makefile
if test 395 -ne `wc -c <Makefile`; then
    echo shar: \"Makefile\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f complex.3 -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"complex.3\"
else
echo shar: Extracting \"complex.3\" \(5830 characters\)
sed "s/^X//" >complex.3 <<'END_OF_complex.3'
X'\" e
X.TH COMPLEX 3V LOCAL
X'\"    last edit:    86/02/03    D A Gwyn
X'\"    SCCS ID:    @(#)complex.3    1.2 (modified for public version)
X.EQ
Xdelim @@
X.EN
X.SH NAME
Xcomplex \- complex arithmetic operations
X.SH SYNOPSIS
X.B
X#include <complex.h>    /* assuming appropriate cc \-I option */
X.br
X/* All the following functions are declared in this header file. */
X.P
X.B complex *CxAdd(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxSub(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxMul(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxDiv(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxSqrt(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxScal(cp,\^s);
X.br
X.B complex *cp;
X.br
X.B double s;
X.P
X.B complex *CxNeg(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxConj(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxCopy(ap,bp);
X.br
X.B complex *ap, *bp;
X.P
X.B complex *CxCons(cp,\^r,\^i);
X.br
X.B complex *cp;
X.br
X.B double r, i;
X.P
X.B complex *CxPhsr(cp,m,p);
X.br
X.B complex *cp;
X.br
X.B double m, p;
X.P
X.B double CxReal(cp);
X.br
X.B complex *cp;
X.P
X.B double CxImag(cp);
X.br
X.B complex *cp;
X.P
X.B double CxAmpl(cp);
X.br
X.B complex *cp;
X.P
X.B double CxPhas(cp);
X.br
X.B complex *cp;
X.P
X.B complex *CxAllo(\ );
X.P
X.B void CxFree(cp);
X.br
X.B complex *cp;
X.SH DESCRIPTION
XThese routines perform arithmetic
Xand other useful operations on complex numbers.
XAn appropriate data structure
X.B complex
Xis defined in the header file;
Xall access to
X.B complex
Xdata should be
X.I via
Xthese predefined functions.
X(See
X.SM HINTS
Xfor further information.)
X.P
XIn the following descriptions,
Xthe names
X.IR a ,
X.IR b ,
Xand
X.I c
Xrepresent the
X.B complex
Xdata addressed by the corresponding pointers
X.IR ap ,
X.IR bp ,
Xand
X.IR cp .
X.P
X.I CxAdd\^
Xadds
X.I b
Xto
X.I a
Xand returns a pointer to the result.
X.P
X.I CxSub
Xsubtracts
X.I b
Xfrom
X.I a
Xand returns a pointer to the result.
X.P
X.I CxMul\^
Xmultiplies
X.I a
Xby
X.I b
Xand returns a pointer to the result.
X.P
X.I CxDiv
Xdivides
X.I a
Xby
X.I b
Xand returns a pointer to the result.
XThe divisor must not be precisely zero.
X.P
X.I CxSqrt
Xreplaces
X.I c
Xby the ``principal value'' of its square root
X(one having a non-negative imaginary part)
Xand returns a pointer to the result.
X.P
X.I CxScal\^
Xmultiplies
X.I c
Xby the scalar
X.I s
Xand returns a pointer to the result.
X.P
X.I CxNeg
Xnegates
X.I c
Xand returns a pointer to the result.
X.P
X.I CxConj
Xconjugates
X.I c
Xand returns a pointer to the result.
X.P
X.I CxCopy
Xassigns the value of
X.I b
Xto
X.I a
Xand returns a pointer to the result.
X.P
X.I CxCons
Xconstructs the complex number
X.I c
Xfrom its real and imaginary parts
X.I r
Xand
X.IR i ,
Xrespectively,
Xand returns a pointer to the result.
X.P
X.I CxPhsr
Xconstructs the complex number
X.I c
Xfrom its ``phasor'' amplitude and phase (given in radians)
X.I m
Xand
X.IR p ,
Xrespectively,
Xand returns a pointer to the result.
X.P
X.I CxReal\^
Xreturns the real part of the complex number
X.IR c .
X.P
X.I CxImag
Xreturns the imaginary part of the complex number
X.IR c .
X.P
X.I CxAmpl\^
Xreturns the amplitude of the complex number
X.IR c .
X.P
X.I CxPhas
Xreturns the phase of the complex number
X.IR c ,
Xas radians in the range @(- pi , pi ]@.
X.P
X.I CxAllo
Xallocates storage for a
X.B complex
Xdatum; it returns
X.SM
X.B NULL
X(defined as 0 in
X.BR <stdio.h> )
Xif not enough storage is available.
X.P
X.I CxFree
Xreleases storage previously allocated by
X.IR CxAllo .
XThe contents of such storage must not be used afterward.
X.SH HINTS
XThe
X.B complex
Xdata type consists of real and imaginary components;
X.I CxReal\^
Xand
X.I CxImag
Xare actually macros that access these components directly.
XThis allows addresses of the components to be taken,
Xas in the following \s-1EXAMPLE\s0.
X.P
XThe complex functions are designed to be nested;
Xsee the following \s-1EXAMPLE\s0.
XFor this reason,
Xmany of them modify the contents of their first parameter.
X.I CxCopy
Xcan be used to create a ``working copy'' of
X.B complex
Xdata that would otherwise be modified.
X.P
XThe square-root function is inherently double-valued;
Xin most applications, both roots should receive equal consideration.
XThe second root is the negative of the ``principal value''.
X.bp
X.SH EXAMPLE
XThe following program is compiled by the command
X.br
X    $ \fIcc \|\-I/usr/local/include \|example.c \|/usr/local/lib/libcomplex.a \|\-lm\fP
X.br
XIt reads in two complex vectors,
Xthen computes and prints their inner product.
X.sp
X.P
X    #include    <stdio.h>
X.br
X    #include    <complex.h>
X.sp
X    main( argc, argv )
X.br
X        int        argc;
X.br
X        char        *argv[\|];
X.br
X        {
X.br
X        int        n;        /* # elements in each array */
X.br
X        int        i;        /* indexes arrays */
X.br
X        complex        a[10], b[10];    /* input vectors */
X.br
X        complex        s;        /* accumulates scalar product */
X.br
X        complex        *c = CxAllo(\|);    /* holds cross-term */
X.sp
X        if ( c == NULL )
X.br
X            {
X.br
X            (void)fprintf( stderr, ``not enough memory\en'' );
X.br
X            return 1;
X.br
X            }
X.br
X        (void)printf( ``\enenter number of elements: '' );
X.br
X        (void)scanf( `` %d'', &n );
X.br
X        /* (There really should be some input validation here.) */
X.br
X        (void) printf( ``\enenter real, imaginary pairs for first array:\en'' );
X.br
X        for ( i = 0; i < n; ++i )
X.br
X            (void)scanf( `` %lg %lg'', &CxReal( &a[i] ), &CxImag( &a[i] ) );
X.br
X        (void)printf( ``\enenter real, imaginary pairs for second array:\en'' );
X.br
X        for ( i = 0; i < n; ++i )
X.br
X            (void)scanf( `` %lg %lg'', &CxReal( &b[i] ), &CxImag( &b[i] ) );
X.br
X        (void)CxCons( &s, 0.0, 0.0 );    /* initialize accumulator */
X.br
X        for ( i = 0; i < n; ++i )
X.br
X            (void)CxAdd( &s, CxMul( &a[i], CxConj( CxCopy( c, &b[i] ) ) ) );
X.br
X        (void)printf( ``\enproduct is (%g,%g)\en'', CxReal( &s ), CxImag( &s ) );
X.br
X        CxFree( c );
X.br
X        return 0;
X.br
X        }
X.SH FILES
X/usr/local/include/complex.h        header file containing definitions
X.br
X/usr/local/lib/libcomplex.a        complex run-time support library
X.SH AUTHORS
XDouglas A. Gwyn, BRL/VLD-VMB
X.br
XJeff Hanes, BRL/VLD-VMB (original version of
X.IR CxSqrt\^ )
END_OF_complex.3
if test 5830 -ne `wc -c <complex.3`; then
    echo shar: \"complex.3\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f complex.h -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"complex.h\"
else
echo shar: Extracting \"complex.h\" \(966 characters\)
sed "s/^X//" >complex.h <<'END_OF_complex.h'
X/*
X    <complex.h> -- definitions for complex arithmetic routines
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)complex.h    1.1 (modified for public version)
X*/
X
X/* "complex number" data type: */
X
Xtypedef struct
X    {
X    double        re;        /* real part */
X    double        im;        /* imaginary part */
X    }    complex;
X
X/* "The future is now": */
X
X#ifdef __STDC__    /* X3J11 */
X#define    _CxGenPtr    void *        /* generic pointer type */
X#else        /* K&R */
X#define    _CxGenPtr    char *        /* generic pointer type */
X#endif
X
X/* functions that are correctly done as macros: */
X
X#define    CxAllo()        ((complex *)malloc( sizeof (complex) ))
X#define    CxFree( cp )        free( (_CxGenPtr)(cp) )
X#define    CxNeg( cp )        CxScal( cp, -1.0 )
X#define    CxReal( cp )        (cp)->re
X#define    CxImag( cp )        (cp)->im
X
Xextern void        free();
Xextern _CxGenPtr    malloc();
X
X/* library functions: */
X
Xextern double    CxAmpl(), CxPhas();
Xextern complex    *CxAdd(), *CxConj(), *CxCons(), *CxCopy(), *CxDiv(),
X        *CxMul(), *CxPhsr(), *CxScal(), *CxSqrt(), *CxSub();
END_OF_complex.h
if test 966 -ne `wc -c <complex.h`; then
    echo shar: \"complex.h\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cx_test.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cx_test.c\"
else
echo shar: Extracting \"cx_test.c\" \(4250 characters\)
sed "s/^X//" >cx_test.c <<'END_OF_cx_test.c'
X/*
X    ctest -- complex arithmetic test
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cx_test.c    1.1 (modified for public version)
X*/
X
X#include    <stdio.h>
X#include    <math.h>
X
X#include    <complex.h>
X
X#define DEGRAD    57.2957795130823208767981548141051703324054724665642
X                    /* degrees per radian */
X#define Abs( x )    ((x) < 0 ? -(x) : (x))
X#define Max( a, b )    ((a) > (b) ? (a) : (b))
X
Xextern void    exit();
X
X#define    Printf    (void)printf
X
X#define    TOL    1.0e-10            /* tolerance for checks */
X
Xstatic int    errs = 0;        /* tally errors */
X
Xstatic void    CCheck(), RCheck();
Xstatic double    RelDif();
X
X
X/*ARGSUSED*/
Xmain( argc, argv )
X    int    argc;
X    char    *argv[];
X    {
X    complex a, *bp, *cp;
X
X    /* CxAllo test */
X    bp = CxAllo();
X    if ( bp == NULL )
X        {
X        Printf( "CxAllo failed\n" );
X        exit( 1 );
X        }
X
X    /* CxReal, CxImag test */
X    CxReal( bp ) = 1.0;
X    CxImag( bp ) = 2.0;
X    RCheck( "CxReal", CxReal( bp ), 1.0 );
X    RCheck( "CxImag", CxImag( bp ), 2.0 );
X
X    /* CxCons test */
X    cp = CxCons( &a, -3.0, -4.0);
X    CCheck( "CxCons 1", a, -3.0, -4.0 );
X    CCheck( "CxCons 2", *cp, -3.0, -4.0 );
X
X    /* CxNeg test */
X    cp = CxNeg( &a );
X    CCheck( "CxNeg 1", a, 3.0, 4.0 );
X    CCheck( "CxNeg 2", *cp, 3.0, 4.0 );
X
X    /* CxCopy test */
X    cp = CxCopy( bp, &a );
X    (void)CxCons( &a, 1.0, sqrt( 3.0 ) );
X    CCheck( "CxCopy 1", *bp, 3.0, 4.0 );
X    CCheck( "CxCopy 2", *cp, 3.0, 4.0 );
X
X    /* CxAmpl, CxPhas test */
X    RCheck( "CxAmpl 1", CxAmpl( &a ), 2.0 );
X    RCheck( "CxPhas 1", CxPhas( &a ) * DEGRAD, 60.0 );
X    /* try other quadrants */
X    a.re = -a.re;
X    RCheck( "CxAmpl 2", CxAmpl( &a ), 2.0 );
X    RCheck( "CxPhas 2", CxPhas( &a ) * DEGRAD, 120.0 );
X    a.im = -a.im;
X    RCheck( "CxAmpl 3", CxAmpl( &a ), 2.0 );
X    RCheck( "CxPhas 3", CxPhas( &a ) * DEGRAD, -120.0 );
X    a.re = -a.re;
X    RCheck( "CxAmpl 4", CxAmpl( &a ), 2.0 );
X    RCheck( "CxPhas 4", CxPhas( &a ) * DEGRAD, -60.0 );
X    /* one more for good measure */
X    RCheck( "CxAmpl 5", CxAmpl( bp ), 5.0 );
X
X    /* CxPhsr test */
X    cp = CxPhsr( &a, 100.0, -20.0 / DEGRAD );
X    RCheck( "CxPhsr 1", CxAmpl( &a ), 100.0 );
X    RCheck( "CxPhsr 2", CxPhas( &a ) * DEGRAD, -20.0 );
X    RCheck( "CxPhsr 3", CxAmpl( cp ), 100.0 );
X    RCheck( "CxPhsr 4", CxPhas( cp ) * DEGRAD, -20.0 );
X
X    /* CxConj test */
X    cp = CxConj( bp );
X    CCheck( "CxConj 1", *bp, 3.0, -4.0 );
X    CCheck( "CxConj 2", *cp, 3.0, -4.0 );
X
X    /* CxScal test */
X    cp = CxScal( bp, 2.0 );
X    CCheck( "CxScal 1", *bp, 6.0, -8.0 );
X    CCheck( "CxScal 2", *cp, 6.0, -8.0 );
X
X    /* CxAdd test */
X    cp = CxAdd( CxCons( &a, -4.0, 11.0 ), bp );
X    CCheck( "CxAdd 1", a, 2.0, 3.0 );
X    CCheck( "CxAdd 2", *cp, 2.0, 3.0 );
X
X    /* CxSub test */
X    cp = CxSub( CxCons( &a, 4.0, 7.0 ), bp );
X    CCheck( "CxSub 1", a, -2.0, 15.0 );
X    CCheck( "CxSub 2", *cp, -2.0, 15.0 );
X
X    /* CxMul test */
X    cp = CxMul( CxCons( bp, -1.0, 3.0 ), CxCons( &a, 1.0, 2.0 ) );
X    CCheck( "CxMul 1", *bp, -7.0, 1.0 );
X    CCheck( "CxMul 2", *cp, -7.0, 1.0 );
X
X    /* CxDiv test */
X    cp = CxDiv( bp, &a );
X    CCheck( "CxDiv 1", *bp, -1.0, 3.0 );
X    CCheck( "CxDiv 2", *cp, -1.0, 3.0 );
X
X    /* CxSqrt and overlapping CxMul tests */
X    (void)CxCons( &a, -1.0, 2.0 );
X    cp = CxSqrt( CxMul( &a, &a ) );
X    CCheck( "CxSqrt 1", a, -1.0, 2.0 );
X    CCheck( "CxSqrt 2", *cp, -1.0, 2.0 );
X    (void)CxCons( &a, 3.0, 2.0 );
X    cp = CxSqrt( CxMul( &a, &a ) );
X    CCheck( "CxSqrt 3", a, 3.0, 2.0 );
X    CCheck( "CxSqrt 4", *cp, 3.0, 2.0 );
X
X    /* CxFree "test" */
X    CxFree( bp );
X
X    return errs;
X    }
X
X
Xstatic void
XCCheck( s, c, r, i )            /* check complex number */
X    char    *s;            /* message string for failure */
X    complex    c;            /* complex to be checked */
X    double    r, i;            /* expected real, imaginary parts */
X    {
X    if ( RelDif( CxReal( &c ), r ) > TOL
X      || RelDif( CxImag( &c ), i ) > TOL
X       )    {
X        ++errs;
X        Printf( "%s; s.b. (%f,%f), was (%g,%g)\n",
X            s, r, i, c.re, c.im
X              );
X        }
X    }
X
X
Xstatic void
XRCheck( s, d, r )            /* check real number */
X    char    *s;            /* message string for failure */
X    double    d;            /* real to be checked */
X    double    r;            /* expected value */
X    {
X    if ( RelDif( d, r ) > TOL )
X        {
X        ++errs;
X        Printf( "%s; s.b. %f, was %g\n", s, r, d );
X        }
X    }
X
X
Xstatic double
XRelDif( a, b )            /* returns relative difference:    */
X    double    a, b;        /* 0.0 if exactly the same,
X                   otherwise ratio of difference
X                   to the larger of the two    */
X    {
X    double    c = Abs( a );
X    double    d = Abs( b );
X
X    d = Max( c, d );
X
X    return d == 0.0 ? 0.0 : Abs( a - b ) / d;
X    }
END_OF_cx_test.c
if test 4250 -ne `wc -c <cx_test.c`; then
    echo shar: \"cx_test.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxadd.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxadd.c\"
else
echo shar: Extracting \"cxadd.c\" \(304 characters\)
sed "s/^X//" >cxadd.c <<'END_OF_cxadd.c'
X/*
X    CxAdd -- add one complex to another
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxadd.c    1.1
X
X    CxAdd( &a, &b )    adds  b  to  a  and returns  &a
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxAdd( ap, bp )
X    register complex    *ap, *bp;    /* may coincide */
X    {
X    ap->re += bp->re;
X    ap->im += bp->im;
X
X    return ap;
X    }
END_OF_cxadd.c
if test 304 -ne `wc -c <cxadd.c`; then
    echo shar: \"cxadd.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxampl.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxampl.c\"
else
echo shar: Extracting \"cxampl.c\" \(278 characters\)
sed "s/^X//" >cxampl.c <<'END_OF_cxampl.c'
X/*
X    CxAmpl -- amplitude (magnitude, modulus, norm) of a complex
X
X    CxAmpl( &c )    returns  |c|
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxampl.c    1.1
X*/
X
X#include    <math.h>
X
X#include    <complex.h>
X
Xdouble
XCxAmpl( cp )
X    register complex    *cp;
X    {
X    return hypot( cp->re, cp->im );
X    }
END_OF_cxampl.c
if test 278 -ne `wc -c <cxampl.c`; then
    echo shar: \"cxampl.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxconj.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxconj.c\"
else
echo shar: Extracting \"cxconj.c\" \(278 characters\)
sed "s/^X//" >cxconj.c <<'END_OF_cxconj.c'
X/*
X    CxConj -- conjugate a complex
X
X    CxConj( &c )    conjugates  c  and returns  &c
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxconj.c    1.1
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxConj( cp )
X    register complex    *cp;
X    {
X    /* (real part unchanged) */
X    cp->im = -cp->im;
X
X    return cp;
X    }
END_OF_cxconj.c
if test 278 -ne `wc -c <cxconj.c`; then
    echo shar: \"cxconj.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxcons.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxcons.c\"
else
echo shar: Extracting \"cxcons.c\" \(329 characters\)
sed "s/^X//" >cxcons.c <<'END_OF_cxcons.c'
X/*
X    CxCons -- construct a complex from real and imaginary parts
X
X    CxCons( &c, re, im )    makes  c = re + i im  and returns  &c
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxcons.c    1.1
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxCons( cp, re, im )
X    register complex    *cp;
X    double            re, im;
X    {
X    cp->re = re;
X    cp->im = im;
X
X    return cp;
X    }
END_OF_cxcons.c
if test 329 -ne `wc -c <cxcons.c`; then
    echo shar: \"cxcons.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxcopy.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxcopy.c\"
else
echo shar: Extracting \"cxcopy.c\" \(264 characters\)
sed "s/^X//" >cxcopy.c <<'END_OF_cxcopy.c'
X/*
X    CxCopy -- copy a complex
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxcopy.c    1.1
X
X    CxCopy( &a, &b )    copies  b  to  a  and returns  &a
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxCopy( ap, bp )
X    complex    *ap, *bp;        /* may coincide */
X    {
X    *ap = *bp;
X
X    return ap;
X    }
END_OF_cxcopy.c
if test 264 -ne `wc -c <cxcopy.c`; then
    echo shar: \"cxcopy.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxdiv.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxdiv.c\"
else
echo shar: Extracting \"cxdiv.c\" \(820 characters\)
sed "s/^X//" >cxdiv.c <<'END_OF_cxdiv.c'
X/*
X    CxDiv -- divide one complex by another
X
X    CxDiv( &a, &b )    divides  a  by  b  and returns  &a;
X            zero divisor fails
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxdiv.c    1.1 (modified for public version)
X*/
X
X#include    <complex.h>
X
X#define Abs( x )    ((x) < 0 ? -(x) : (x))
X
Xcomplex *
XCxDiv( ap, bp )
X    register complex    *ap, *bp;    /* may coincide (?) */
X    {
X    double    r, s;
X    double    ap__re = ap->re;
X
X    /* Note: classical formula may cause unnecessary overflow */
X    r = bp->re;
X    s = bp->im;
X    if ( Abs( r ) >= Abs( s ) )
X        {
X        r = s / r;        /* <= 1 */
X        s = bp->re + r * s;
X        ap->re = (ap->re + ap->im * r) / s;
X        ap->im = (ap->im - ap__re * r) / s;
X        }
X    else /* Abs( s ) > Abs( r ) */
X        {
X        r = r / s;        /* < 1 */
X        s = s + r * bp->re;
X        ap->re = (ap->re * r + ap->im) / s;
X        ap->im = (ap->im * r - ap__re) / s;
X        }
X
X    return ap;
X    }
END_OF_cxdiv.c
if test 820 -ne `wc -c <cxdiv.c`; then
    echo shar: \"cxdiv.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxmul.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxmul.c\"
else
echo shar: Extracting \"cxmul.c\" \(424 characters\)
sed "s/^X//" >cxmul.c <<'END_OF_cxmul.c'
X/*
X    CxMul -- multiply one complex by another
X
X    CxMul( &a, &b )    multiplies  a  by  b  and returns  &a
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxmul.c    1.1
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxMul( ap, bp )
X    register complex    *ap, *bp;    /* (may coincide) */
X    {
X    double            ap__re = ap->re;
X    double            bp__re = bp->re;
X
X    ap->re = ap__re * bp__re - ap->im * bp->im;
X    ap->im = ap__re * bp->im + ap->im * bp__re;
X
X    return ap;
X    }
END_OF_cxmul.c
if test 424 -ne `wc -c <cxmul.c`; then
    echo shar: \"cxmul.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxphas.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxphas.c\"
else
echo shar: Extracting \"cxphas.c\" \(406 characters\)
sed "s/^X//" >cxphas.c <<'END_OF_cxphas.c'
X/*
X    CxPhas -- phase (angle, argument) of a complex
X
X    CxPhas( &c )    returns  arg(c)  in radians (-Pi,Pi]
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxphas.c    1.1 (modified for public version)
X*/
X
X#include    <math.h>
X
X#include    <complex.h>
X
Xdouble
XCxPhas( cp )
X    register complex    *cp;
X    {
X    if ( cp->re == 0.0 && cp->im == 0.0 )
X        return 0.0;        /* can't trust atan2() */
X    else
X        return atan2( cp->im, cp->re );
X    }
END_OF_cxphas.c
if test 406 -ne `wc -c <cxphas.c`; then
    echo shar: \"cxphas.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxphsr.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxphsr.c\"
else
echo shar: Extracting \"cxphsr.c\" \(395 characters\)
sed "s/^X//" >cxphsr.c <<'END_OF_cxphsr.c'
X/*
X    CxPhsr -- construct a complex "phasor" from amplitude and phase
X
X    CxPhsr( &c, amp, phs )    makes  c = amp exp(i phs)
X                    and returns  &c
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxphsr.c    1.1
X*/
X
X#include    <math.h>
X
X#include    <complex.h>
X
Xcomplex *
XCxPhsr( cp, amp, phs )
X    register complex    *cp;
X    double            amp, phs;
X    {
X    cp->re = amp * cos( phs );
X    cp->im = amp * sin( phs );
X
X    return cp;
X    }
END_OF_cxphsr.c
if test 395 -ne `wc -c <cxphsr.c`; then
    echo shar: \"cxphsr.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxscal.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxscal.c\"
else
echo shar: Extracting \"cxscal.c\" \(291 characters\)
sed "s/^X//" >cxscal.c <<'END_OF_cxscal.c'
X/*
X    CxScal -- multiply a complex by a scalar
X
X    CxScal( &c, s )    scales  c  by  s  and returns  &c
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxscal.c    1.1
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxScal( cp, s )
X    register complex    *cp;
X    double            s;
X    {
X    cp->re *= s;
X    cp->im *= s;
X
X    return cp;
X    }
END_OF_cxscal.c
if test 291 -ne `wc -c <cxscal.c`; then
    echo shar: \"cxscal.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxsqrt.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxsqrt.c\"
else
echo shar: Extracting \"cxsqrt.c\" \(1339 characters\)
sed "s/^X//" >cxsqrt.c <<'END_OF_cxsqrt.c'
X/*
X    CxSqrt -- compute square root of complex number
X
X    CxSqrt( &c )    replaces  c  by  sqrt(c)  and returns  &c
X
X    Note:    This is a double-valued function; the result of
X        CxSqrt() always has nonnegative imaginary part.
X
X    inspired by Jeff Hanes' version
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxsqrt.c    1.1 (modified for public version)
X*/
X
X#include    <math.h>
X
X#include    <complex.h>
X
X#define    Sgn( x )    ((x) == 0 ? 0 : (x) > 0 ? 1 : -1)
X
Xcomplex *
XCxSqrt( cp )
X    register complex    *cp;
X    {
X    /* record signs of original real & imaginary parts */
X    int            re_sign = Sgn( cp->re );
X    int            im_sign = Sgn( cp->im );
X
X    /* special cases are not necessary; they are here for speed */
X
X    if ( re_sign == 0 )
X        if ( im_sign == 0 )
X            ;        /* (0,0) already there */
X        else if ( im_sign > 0 )
X            cp->re = cp->im = sqrt( cp->im / 2.0 );
X        else            /* im_sign < 0 */
X            cp->re = -(cp->im = sqrt( -cp->im / 2.0 ));
X    else if ( im_sign == 0 )
X        if ( re_sign > 0 )
X            cp->re = sqrt( cp->re );
X/*            cp->im = 0.0;    /* 0 already there */
X        else    {        /* re_sign < 0 */
X            cp->im = sqrt( -cp->re );
X            cp->re = 0.0;
X            }
X    else    {            /* no shortcuts */
X        double    ampl = CxAmpl( cp );
X
X        cp->im = sqrt( (ampl - cp->re) /2.0 );
X
X        if ( im_sign > 0 )
X            cp->re = sqrt( (ampl + cp->re) / 2.0 );
X        else            /* im_sign < 0 */
X            cp->re = -sqrt( (ampl + cp->re) / 2.0 );
X        }
X
X    return cp;
X    }
END_OF_cxsqrt.c
if test 1339 -ne `wc -c <cxsqrt.c`; then
    echo shar: \"cxsqrt.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f cxsub.c -a "${1}" != "-c" ; then
  echo shar: Will not over-write existing file \"cxsub.c\"
else
echo shar: Extracting \"cxsub.c\" \(320 characters\)
sed "s/^X//" >cxsub.c <<'END_OF_cxsub.c'
X/*
X    CxSub -- subtract one complex from another
X
X    CxSub( &a, &b )    subtracts  b  from  a  and returns  &a
X
X    last edit:    86/01/04    D A Gwyn
X
X    SCCS ID:    @(#)cxsub.c    1.1
X*/
X
X#include    <complex.h>
X
Xcomplex *
XCxSub( ap, bp )
X    register complex    *ap, *bp;    /* (may coincide) */
X    {
X    ap->re -= bp->re;
X    ap->im -= bp->im;
X
X    return ap;
X    }
END_OF_cxsub.c
if test 320 -ne `wc -c <cxsub.c`; then
    echo shar: \"cxsub.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of shell archive.
exit 0
--

Rich $alz            "Anger is an energy"
Cronus Project, BBN Labs    [e-mail]
Moderator, comp.sources.unix    [e-mail]
--

Rich $alz            "Anger is an energy"
Cronus Project, BBN Labs    [e-mail]
Moderator, comp.sources.unix    [e-mail]